package fem2.material;

import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

import java.util.ArrayList;

import math2.MathUtilities;
import math2.MutableDouble;
import fem2.Debugger;
import fem2.Element;
import fem2.Variable;
import fem2.element.StructuralElement;
import fem2.enu.EchoLevelType;
import fem2.enu.State;
import fem2.enu.VariableName;

/**
 * Isotropic material model using energy norm as equivalent strain measure. For
 * detail please refer to
 * 
 * <pre>
 * An implicit/explicit integration scheme to increase computabilty of nonlinear material and contact/friction problem
 * </pre>
 * 
 * by J. Oliver
 * 
 * @author hbui
 * 
 */
public class IsotropicDamageMM extends InelasticMaterialModel {

	protected FluencyCriteria softeningRule;
	protected double alpha0;

	/**
	 * 
	 * @param E
	 *            elastic modulus
	 * @param nu
	 *            Poisson ratio
	 * @param rho
	 *            density
	 * @param thickness
	 *            thickness
	 * @param ss
	 *            stress state
	 * @param softeningRule
	 *            softening rule
	 * @param alpha0
	 *            the initial value of internal variable and stress-like
	 *            variable for this model. alpha0 must satisfy
	 *            softeningRule(alpha0) = alpha0 to ensure damage is 0 in the
	 *            beginning. Normally, alpha0 = ft/sqrt(E)
	 */
	public IsotropicDamageMM(double E, double nu, double rho, double thickness, State ss,
			FluencyCriteria softeningRule, double alpha0) {
		super(E, nu, thickness, rho, ss);

		this.softeningRule = softeningRule;
		this.alpha0 = alpha0;

		/*
		 * check if initial softening condition is satisfied
		 */
		if (Math.abs(softeningRule.computeFluencyValue(alpha0) - alpha0) > MathUtilities.tol12) {
			throw new Error("initial softening condition is not satisfied");
		}
	}

	@Override
	public void Initialize(Element e) {
		super.Initialize(e);

		ArrayList<Variable> vars = getVariables(e);
		for (Variable v : vars) {
			v.registerScalarVariable(VariableName.ALPHA, true);
			v.registerScalarVariable(VariableName.STRESS_LIKE_VARIABLE, true);// q
			v.registerScalarVariable(VariableName.DAMAGE, false);
			/*
			 * initialize internal variable and stress-like variable
			 */
			MutableDouble alpha = v.getScalarVariable(VariableName.ALPHA);
			MutableDouble q = v.getScalarVariable(VariableName.STRESS_LIKE_VARIABLE);
			alpha.setValue(alpha0);
			q.setValue(alpha0);
			/*
			 * save the variable to create history
			 */
			v.save();
		}
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
		StructuralElement se = (StructuralElement) e;
		int ng = e.getNumIntegrationPoints();

		for (int k = 0; k < ng; k++) {
			Variable var = getVariables(e).get(k);
			double[] xi = e.getIntegrationPoint(k);
			IVector sigma = se.getStress(k);
			IMatrix C = getConstitutiveMatrix(e, k);
			IVector epsilon = var.getVectorVariable(VariableName.STRAIN);
			MutableDouble alpha = var.getScalarVariable(VariableName.ALPHA);
			double alpha_old = var.getScalarVariable(VariableName.ALPHA, -1);
			MutableDouble q = var.getScalarVariable(VariableName.STRESS_LIKE_VARIABLE);
			double q_old = var.getScalarVariable(VariableName.STRESS_LIKE_VARIABLE, -1);
			MutableDouble damage = var.getScalarVariable(VariableName.DAMAGE);

			/*
			 * compute strain
			 */
			BLAM.copy(se.computeStrain(xi), epsilon);

			/*
			 * compute effective stress
			 */
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, De, epsilon, 0.0, sigma);

			/*
			 * compute equivalent strain
			 */
			double etilde = Math.sqrt(BLAM.dot(epsilon, sigma));

			if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
				Debugger.watch("etilde = ", etilde);
			}

			/*
			 * compute trial state
			 */
			double alpha_trial = alpha_old;

			/*
			 * check loading/unloading criteria (use g to be consistent with the
			 * paper)
			 */
			double g = etilde - alpha_trial;
			double dlambda;
			if (g <= 0.0) {
				/*
				 * unloading/elastic state
				 */
				dlambda = 0.0;
			} else {
				/*
				 * loading/damage state
				 */
				dlambda = g;
			}

			/*
			 * update internal variable
			 */
			alpha.setValue(alpha_old + dlambda);

			/*
			 * update stress-like variable
			 */
			double H = softeningRule.computeFluencyDerivative(alpha_old);
			// double H =
			// softeningRule.computeFluencyDerivative(alpha.getValue());
			q.setValue(q_old + H * dlambda);

			/*
			 * update damage variable
			 */
			damage.setValue(1 - q.getValue() / alpha.getValue());
			// if (damage.getValue() > 0.9999) {
			// damage.setValue(0.9999);
			// }
			if (damage.getValue() < 0.0) {
				/*
				 * just for debugging, to make sure the input & material model
				 * work correcly
				 */
				throw new Error("negative damage. It's unphysical");
			}

			/*
			 * update tangent operator
			 */
			if (dlambda <= 0.0) {
				BLAM.copy(De, C);
				BLAM.scale(1 - damage.getValue(), C);
			} else {
				BLAM.copy(De, C);
				BLAM.scale(1 - damage.getValue(), C);
				BLAM.multiply(
						-(q.getValue() - H * alpha.getValue()) / Math.pow(alpha.getValue(), 3),
						sigma, sigma, C);
			}

			/*
			 * update stress
			 */
			BLAM.scale(1 - damage.getValue(), sigma);

			if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
				Debugger.watch("finalize step:");
				Debugger.watch("epsilon =" + epsilon);
				Debugger.watch("sigma =" + sigma);
				Debugger.watch("alpha = ", var.getScalarVariable(VariableName.ALPHA));
				Debugger.watch("damage = ", var.getScalarVariable(VariableName.DAMAGE));
			}
		}
	}

}
